library(rgrass7)
## Loading required package: XML
## GRASS GIS interface loaded with GRASS version: GRASS 7.8.2 (2019)
## and location: pantuflas
gisdbase <- 'grass-data-test' #Base de datos de GRASS GIS
wd <- getwd() #Directorio de trabajo
wd
## [1] "/home/jr/unidad-4-asignacion-1-procesos-fluviales"
loc <- initGRASS(gisBase = "/usr/lib/grass78/",
                 home = wd,
                 gisDbase = paste(wd, gisdbase, sep = '/'),
                 location = 'pantuflas',
                 mapset = "PERMANENT",
                 override = TRUE)
gmeta()
## gisdbase    /home/jr/unidad-4-asignacion-1-procesos-fluviales/grass-data-test 
## location    pantuflas 
## mapset      PERMANENT 
## rows        892 
## columns     971 
## north       2122920 
## south       2042640 
## west        253800 
## east        341190 
## nsres       90 
## ewres       90 
## projection  +proj=utm +no_defs +zone=19 +a=6378137 +rf=298.257223563
## +towgs84=0.000,0.000,0.000 +type=crs +to_meter=1
execGRASS(
  'g.list',
  flags = 't',
  parameters = list(
    type = c('raster', 'vector')
  )
)
## raster/demint
## raster/rbasin_pant_demint_accumulation
## raster/rbasin_pant_demint_aspect
## raster/rbasin_pant_demint_dist2out
## raster/rbasin_pant_demint_drainage
## raster/rbasin_pant_demint_hack
## raster/rbasin_pant_demint_hillslope_distance
## raster/rbasin_pant_demint_horton
## raster/rbasin_pant_demint_shreve
## raster/rbasin_pant_demint_slope
## raster/rbasin_pant_demint_strahler
## raster/stream-de-rstr
## vector/rbasin_pant_demint_basin
## vector/rbasin_pant_demint_mainchannel
## vector/rbasin_pant_demint_network
## vector/rbasin_pant_demint_outlet
## vector/rbasin_pant_demint_outlet_snap
## vector/stream_de_rstr

Convertir a números enteros la extensión y la resolución del DEM

library(raster)
## Loading required package: sp
rutadem <- 'data/dem.tif'
rawextent <- extent(raster(rutadem))
rawextent
## class      : Extent 
## xmin       : 253740.1 
## xmax       : 341262.4 
## ymin       : 2042609 
## ymax       : 2122949
devtools::source_url('https://raw.githubusercontent.com/geofis/rgrass/master/integerextent.R')
## SHA-1 hash of file is 65b57fe7cb20cd1976572cd0a7e98def9e95d83c
devtools::source_url('https://raw.githubusercontent.com/geofis/rgrass/master/xyvector.R')
## SHA-1 hash of file is 89fa5ae436d9a7d7a0c799b789c560eb5e421cfd
newextent <- intext(e = rawextent, r = 90, type = 'inner')
newextent
## class      : Extent 
## xmin       : 253800 
## xmax       : 341190 
## ymin       : 2042640 
## ymax       : 2122920
gdalUtils::gdalwarp(
  srcfile = 'data/dem.tif',
  dstfile = 'data/demint.tif',
  te = xyvector(newextent),
  tr = c(90,90),
  r = 'bilinear',
  overwrite = T
)
## Registered S3 method overwritten by 'R.oo':
##   method        from       
##   throw.default R.methodsS3
## NULL

Importar a sesión de GRASS

rutademint <- 'data/demint.tif'
execGRASS(
  "g.proj",
  flags = c('t','c'),
  georef=rutademint)
gmeta()
execGRASS(
  "r.in.gdal",
  flags='overwrite',
  parameters=list(
    input=rutademint,
    output="demint"
  )
)
execGRASS(
  "g.region",
  parameters=list(
    raster = "demint",
    align = "demint"
  )
)
gmeta()
## gisdbase    /home/jr/unidad-4-asignacion-1-procesos-fluviales/grass-data-test 
## location    pantuflas 
## mapset      PERMANENT 
## rows        892 
## columns     971 
## north       2122920 
## south       2042640 
## west        253800 
## east        341190 
## nsres       90 
## ewres       90 
## projection  +proj=utm +no_defs +zone=19 +a=6378137 +rf=298.257223563
## +towgs84=0.000,0.000,0.000 +type=crs +to_meter=1
execGRASS(
  'g.list',
  flags = 't',
  parameters = list(
    type = c('raster', 'vector')
  )
)
## raster/demint
## raster/rbasin_pant_demint_accumulation
## raster/rbasin_pant_demint_aspect
## raster/rbasin_pant_demint_dist2out
## raster/rbasin_pant_demint_drainage
## raster/rbasin_pant_demint_hack
## raster/rbasin_pant_demint_hillslope_distance
## raster/rbasin_pant_demint_horton
## raster/rbasin_pant_demint_shreve
## raster/rbasin_pant_demint_slope
## raster/rbasin_pant_demint_strahler
## raster/stream-de-rstr
## vector/rbasin_pant_demint_basin
## vector/rbasin_pant_demint_mainchannel
## vector/rbasin_pant_demint_network
## vector/rbasin_pant_demint_outlet
## vector/rbasin_pant_demint_outlet_snap
## vector/stream_de_rstr

Generar red de drenaje para obtener coordenada posteriormente

execGRASS(
  "r.stream.extract",
  flags = c('overwrite','quiet'),
  parameters = list(
    elevation = 'demint',
    threshold = 80,
    stream_raster = 'stream-de-rstr',
    stream_vector = 'stream_de_rstr'
  )
)
## Warning in execGRASS("r.stream.extract", flags = c("overwrite", "quiet"), : The command:
## r.stream.extract --overwrite --quiet elevation=demint threshold=80 stream_raster=stream-de-rstr stream_vector=stream_de_rstr
## produced at least one warning during execution:
## WARNING: Vector map <stream_de_rstr> already exists and will be overwritten
## WARNING: Vector map <stream_de_rstr> already exists and will be overwritten
execGRASS(
  'g.list',
  flags = 't',
  parameters = list(
    type = c('raster', 'vector')
  )
)
## raster/demint
## raster/rbasin_pant_demint_accumulation
## raster/rbasin_pant_demint_aspect
## raster/rbasin_pant_demint_dist2out
## raster/rbasin_pant_demint_drainage
## raster/rbasin_pant_demint_hack
## raster/rbasin_pant_demint_hillslope_distance
## raster/rbasin_pant_demint_horton
## raster/rbasin_pant_demint_shreve
## raster/rbasin_pant_demint_slope
## raster/rbasin_pant_demint_strahler
## raster/stream-de-rstr
## vector/rbasin_pant_demint_basin
## vector/rbasin_pant_demint_mainchannel
## vector/rbasin_pant_demint_network
## vector/rbasin_pant_demint_outlet
## vector/rbasin_pant_demint_outlet_snap
## vector/stream_de_rstr

Obtener coordenada

library(sp)
use_sp()
library(mapview)
netw <- spTransform(
  readVECT('stream_de_rstr'),
  CRSobj = CRS("+init=epsg:4326"))
mapview(netw, col.regions = 'blue', legend = FALSE)
## Warning in CPL_crs_from_input(x): GDAL Message 1: +init=epsg:XXXX syntax is
## deprecated. It might return a CRS with a non-EPSG compliant axis order.
## Warning: sf layer has inconsistent datum (+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs).
## Need '+proj=longlat +datum=WGS84'

Transformar coordenada a EPSG:32619 como número entero

source('my-trans.R')
outlet <- as.integer(my_trans(c(-70.77398,18.90123)))

Ejecutar r.basin

pref <- 'rbasin_pant'
execGRASS(
  "r.basin",
  flags = 'overwrite',
  parameters = list(
    map = 'demint',
    prefix = pref,
    coordinates = outlet,
    threshold = 80,
    dir = 'salidas-rbasin/pantuflas'
  )
)
execGRASS(
  'g.list',
  flags = 't',
  parameters = list(
    type = c('raster', 'vector')
  )
)
## raster/demint
## raster/rbasin_pant_demint_accumulation
## raster/rbasin_pant_demint_aspect
## raster/rbasin_pant_demint_dist2out
## raster/rbasin_pant_demint_drainage
## raster/rbasin_pant_demint_hack
## raster/rbasin_pant_demint_hillslope_distance
## raster/rbasin_pant_demint_horton
## raster/rbasin_pant_demint_shreve
## raster/rbasin_pant_demint_slope
## raster/rbasin_pant_demint_strahler
## raster/stream-de-rstr
## vector/rbasin_pant_demint_basin
## vector/rbasin_pant_demint_mainchannel
## vector/rbasin_pant_demint_network
## vector/rbasin_pant_demint_outlet
## vector/rbasin_pant_demint_outlet_snap
## vector/stream_de_rstr

Si r.basin arrojara error (sólo en el caso de error, no en caso de advertencia), ejecutar este bloque para borrar las salidas anteriores y reejecutar el r.basin:

execGRASS(
  "g.remove",
  flags = 'f',
  parameters = list(
    type = c('raster','vector'),
    pattern = paste0(pref, '*')
  )
)

Cargar los vectoriales transformados a EPSG:4326 para visualizar en leaflet

rbnetw <- spTransform(
  readVECT('rbasin_pant_demint_network'),
  CRSobj = CRS("+init=epsg:4326"))
rbnetw
rbmain <- spTransform(
  readVECT('rbasin_pant_demint_mainchannel'),
  CRSobj = CRS("+init=epsg:4326"))
rbmain
rbbasin <- spTransform(
  readVECT('rbasin_pant_demint_basin'),
  CRSobj = CRS("+init=epsg:4326"))
rbbasin
library(leaflet)
leaflet() %>%
  addProviderTiles(providers$Stamen.Terrain, group = 'terrain') %>%
  addPolylines(data = rbnetw, weight = 3, opacity = 0.7) %>% 
  addPolylines(data = rbmain, weight = 3, opacity = 0.7, color = 'red') %>% 
  addPolygons(data = rbbasin) %>% 
  leafem::addHomeButton(extent(rbbasin), 'Ver cuenca')

Explorar los parámetros de cuenca

library(readr)
rbpantpar1 <- read_csv("salidas-rbasin/pantuflas/rbasin_pant_demint_parametersT.csv")
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   Rectangle_containing_basin_N_W = col_character(),
##   Rectangle_containing_basin_S_E = col_character()
## )
## See spec(...) for full column specifications.
rbpantpar1 %>% tibble::as_tibble()
## # A tibble: 1 x 34
##        x      y Easting_Centroi… Northing_Centro… Rectangle_conta…
##    <dbl>  <dbl>            <dbl>            <dbl> <chr>           
## 1 313155 2.09e6           319815          2087685 ('313110', '209…
## # … with 29 more variables: Rectangle_containing_basin_S_E <chr>,
## #   Area_of_basin_km2 <dbl>, Perimeter_of_basin_km <dbl>,
## #   Max_Elevation <dbl>, Min_Elevation <dbl>, Elevation_Difference <dbl>,
## #   Mean_Elevation <dbl>, Mean_Slope <dbl>,
## #   Length_of_Directing_Vector_km <dbl>,
## #   Prevalent_Orientation_deg_from_north_ccw <dbl>,
## #   Compactness_Coefficient <dbl>, Circularity_Ratio <dbl>,
## #   Topological_Diameter <dbl>, Elongation_Ratio <dbl>,
## #   Shape_Factor <dbl>, Concentration_Time_hr <dbl>,
## #   Length_of_Mainchannel_km <dbl>,
## #   Mean_slope_of_mainchannel_percent <dbl>,
## #   Mean_hillslope_length_m <dbl>, Magnitudo <dbl>,
## #   Max_order_Strahler <dbl>, Number_of_streams <dbl>,
## #   Total_Stream_Length_km <dbl>, First_order_stream_frequency <dbl>,
## #   Drainage_Density_km_over_km2 <dbl>, Bifurcation_Ratio_Horton <dbl>,
## #   Length_Ratio_Horton <dbl>, Area_ratio_Horton <dbl>,
## #   Slope_ratio_Horton <dbl>

rbpantpar2 <- read_csv(
  "salidas-rbasin/pantuflas/rbasin_pant_demint_parameters.csv",
  skip=2, col_names = c('Parameter', 'Value'))
## Parsed with column specification:
## cols(
##   Parameter = col_character(),
##   Value = col_character()
## )
rbpantpar2 %>% print(n=Inf)
## # A tibble: 32 x 2
##    Parameter                                             Value             
##    <chr>                                                 <chr>             
##  1 Easting Centroid of basin                             319815.00         
##  2 Northing Centroid of basin                            2087685.00        
##  3 Rectangle containing basin N-W                        ('313110', '20975…
##  4 Rectangle containing basin S-E                        ('326880', '20782…
##  5 Area of basin [km^2]                                  128.6087625       
##  6 Perimeter of basin [km]                               63.9170486606335  
##  7 Max Elevation [m s.l.m.]                              2619.45           
##  8 Min Elevation [m s.l.m.]                              1054.484          
##  9 Elevation Difference [m]                              1564.966          
## 10 Mean Elevation                                        1588.877          
## 11 Mean Slope                                            13.74             
## 12 Length of Directing Vector [km]                       7.367794853278693 
## 13 Prevalent Orientation [degree from north, counterclo… 0.441915859817346…
## 14 Compactness Coefficient                               4.994895129597467 
## 15 Circularity Ratio                                     0.395591541103538…
## 16 Topological Diameter                                  18.0              
## 17 Elongation Ratio                                      0.6185502219049736
## 18 Shape Factor                                          6.216632397877906 
## 19 Concentration Time (Giandotti, 1934) [hr]             2.413889375890686 
## 20 Length of Mainchannel [km]                            20.687850635      
## 21 Mean slope of mainchannel [percent]                   6.914012642018835 
## 22 Mean hillslope length [m]                             858.9935          
## 23 Magnitudo                                             33.0              
## 24 Max order (Strahler)                                  4                 
## 25 Number of streams                                     44                
## 26 Total Stream Length [km]                              99.8046           
## 27 First order stream frequency                          0.2565921587185787
## 28 Drainage Density [km/km^2]                            0.7760326595164928
## 29 Bifurcation Ratio (Horton)                            3.2429            
## 30 Length Ratio (Horton)                                 0.8840            
## 31 Area ratio (Horton)                                   3.7251            
## 32 Slope ratio (Horton)                                  1.5582

Limpiar archivo de bloqueo del conjunto de mapas de GRASS

unlink_.gislock()

Referencias